home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-04 / amsf20.zip / EIGEN.FOR < prev    next >
Text File  |  1992-01-06  |  5KB  |  196 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE EIGEN
  5. C
  6. C        PURPOSE
  7. C           COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
  8. C           MATRIX
  9. C
  10. C        USAGE
  11. C           CALL EIGEN(A,R,N,MV)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION.
  15. C               RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF
  16. C               MATRIX A IN DESCENDING ORDER.
  17. C           R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE,
  18. C               IN SAME SEQUENCE AS EIGENVALUES)
  19. C           N - ORDER OF MATRICES A AND R
  20. C           MV- INPUT CODE
  21. C                   0   COMPUTE EIGENVALUES AND EIGENVECTORS
  22. C                   1   COMPUTE EIGENVALUES ONLY (R NEED NOT BE
  23. C                       DIMENSIONED BUT MUST STILL APPEAR IN CALLING
  24. C                       SEQUENCE)
  25. C
  26. C        REMARKS
  27. C           ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1)
  28. C           MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R
  29. C
  30. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  31. C           NONE
  32. C
  33. C        METHOD
  34. C           DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED
  35. C           BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL
  36. C           METHODS FOR DIGITAL COMPUTERS', EDITED BY A. RALSTON AND
  37. C           H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7
  38. C
  39. C     ..................................................................
  40. C
  41.       SUBROUTINE EIGEN(A,R,N,MV)
  42.       DIMENSION A(1),R(1)
  43. C
  44. C        ...............................................................
  45. C
  46. C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
  47. C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
  48. C        STATEMENT WHICH FOLLOWS.
  49. C
  50.       DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,
  51.      1                 COSX2,SINCS,RANGE
  52. C
  53. C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
  54. C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
  55. C        ROUTINE.
  56. C
  57. C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
  58. C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENTS
  59. C        40, 68, 75, AND 78 MUST BE CHANGED TO DSQRT.  ABS IN STATEMENT
  60. C        62 MUST BE CHANGED TO DABS. THE CONSTANT IN STATEMENT 5 SHOULD
  61. C        BE CHANGED TO 1.0D-12.
  62. C
  63. C        ...............................................................
  64. C
  65. C        GENERATE IDENTITY MATRIX
  66. C
  67.     5 RANGE=1.0E-6
  68.       IF(MV-1) 10,25,10
  69.    10 IQ=-N
  70.       DO 20 J=1,N
  71.       IQ=IQ+N
  72.       DO 20 I=1,N
  73.       IJ=IQ+I
  74.       R(IJ)=0.0
  75.       IF(I-J) 20,15,20
  76.    15 R(IJ)=1.0
  77.    20 CONTINUE
  78. C
  79. C        COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX)
  80. C
  81.    25 ANORM=0.0
  82.       DO 35 I=1,N
  83.       DO 35 J=I,N
  84.       IF(I-J) 30,35,30
  85.    30 IA=I+(J*J-J)/2
  86.       ANORM=ANORM+A(IA)*A(IA)
  87.    35 CONTINUE
  88.       IF(ANORM) 165,165,40
  89.    40 ANORM=1.414*SQRT(ANORM)
  90.       ANRMX=ANORM*RANGE/FLOAT(N)
  91. C
  92. C        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
  93. C
  94.       IND=0
  95.       THR=ANORM
  96.    45 THR=THR/FLOAT(N)
  97.    50 L=1
  98.    55 M=L+1
  99. C
  100. C        COMPUTE SIN AND COS
  101. C
  102.    60 MQ=(M*M-M)/2
  103.       LQ=(L*L-L)/2
  104.       LM=L+MQ
  105.    62 IF( ABS(A(LM))-THR) 130,65,65
  106.    65 IND=1
  107.       LL=L+LQ
  108.       MM=M+MQ
  109.       X=0.5*(A(LL)-A(MM))
  110.    68 Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X)
  111.       IF(X) 70,75,75
  112.    70 Y=-Y
  113.    75 SINX=Y/ SQRT(2.0*(1.0+( SQRT(1.0-Y*Y))))
  114.       SINX2=SINX*SINX
  115.    78 COSX= SQRT(1.0-SINX2)
  116.       COSX2=COSX*COSX
  117.       SINCS =SINX*COSX
  118. C
  119. C        ROTATE L AND M COLUMNS
  120. C
  121.       ILQ=N*(L-1)
  122.       IMQ=N*(M-1)
  123.       DO 125 I=1,N
  124.       IQ=(I*I-I)/2
  125.       IF(I-L) 80,115,80
  126.    80 IF(I-M) 85,115,90
  127.    85 IM=I+MQ
  128.       GO TO 95
  129.    90 IM=M+IQ
  130.    95 IF(I-L) 100,105,105
  131.   100 IL=I+LQ
  132.       GO TO 110
  133.   105 IL=L+IQ
  134.   110 X=A(IL)*COSX-A(IM)*SINX
  135.       A(IM)=A(IL)*SINX+A(IM)*COSX
  136.       A(IL)=X
  137.   115 IF(MV-1) 120,125,120
  138.   120 ILR=ILQ+I
  139.       IMR=IMQ+I
  140.       X=R(ILR)*COSX-R(IMR)*SINX
  141.       R(IMR)=R(ILR)*SINX+R(IMR)*COSX
  142.       R(ILR)=X
  143.   125 CONTINUE
  144.       X=2.0*A(LM)*SINCS
  145.       Y=A(LL)*COSX2+A(MM)*SINX2-X
  146.       X=A(LL)*SINX2+A(MM)*COSX2+X
  147.       A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
  148.       A(LL)=Y
  149.       A(MM)=X
  150. C
  151. C        TESTS FOR COMPLETION
  152. C
  153. C        TEST FOR M = LAST COLUMN
  154. C
  155.   130 IF(M-N) 135,140,135
  156.   135 M=M+1
  157.       GO TO 60
  158. C
  159. C        TEST FOR L = SECOND FROM LAST COLUMN
  160. C
  161.   140 IF(L-(N-1)) 145,150,145
  162.   145 L=L+1
  163.       GO TO 55
  164.   150 IF(IND-1) 160,155,160
  165.   155 IND=0
  166.       GO TO 50
  167. C
  168. C        COMPARE THRESHOLD WITH FINAL NORM
  169. C
  170.   160 IF(THR-ANRMX) 165,165,45
  171. C
  172. C        SORT EIGENVALUES AND EIGENVECTORS
  173. C
  174.   165 IQ=-N
  175.       DO 185 I=1,N
  176.       IQ=IQ+N
  177.       LL=I+(I*I-I)/2
  178.       JQ=N*(I-2)
  179.       DO 185 J=I,N
  180.       JQ=JQ+N
  181.       MM=J+(J*J-J)/2
  182.       IF(A(LL)-A(MM)) 170,185,185
  183.   170 X=A(LL)
  184.       A(LL)=A(MM)
  185.       A(MM)=X
  186.       IF(MV-1) 175,185,175
  187.   175 DO 180 K=1,N
  188.       ILR=IQ+K
  189.       IMR=JQ+K
  190.       X=R(ILR)
  191.       R(ILR)=R(IMR)
  192.   180 R(IMR)=X
  193.   185 CONTINUE
  194.       RETURN
  195.       END
  196.